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The Hubbard model with additional intersite interaction l V ( the extended Hubbard 
model) is investigated by the correlator projection method (CPM). CPM is a newly devel- 
oped numerical method which combines the equation-of-motion approach and the dynamical 
mean-field theory. Using this method, properties of the extended Hubbard Model at quar- 
ter filling are discussed with special emphasis on the metal-insulator transition induced by 
electron-electron correlations. The result shows a metal-insulator transition to a charge or- 
dered insulator with antiferromagnetic order at low temperatures, but a charge ordered insu- 
lator without magnetic symmetry breaking at intermediate temperatures. Here, the magnetic 
order is found to be suppressed to low temperatures because of the smallness of the exchange 
coupling J c ff. The present results are in sharp contrast with the Hatree-Fock approximation 
whereas in agreement with the experimental results on quarter-filled materials with strong 
correlations, such as organic conductors. 

KEYWORDS: extended Hubbard model, metal-insulator transition, intersite Coulomb interac- 
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1. Introduction 

Translational symmetry breaking in the charge degrees of freedom, or charge ordering, is 
widely observed in strongly correlated electron systems such as manganites, 1-3 vanadates, 4 and 
various organic conductors. 5 ' 6 In these systems, the charge-ordering transition often accom- 
pany a large jump in electric resistivity or metal-insulator transition (MIT). 7 It has attracted 
much interest both from theoretical and experimental point of views. 7-11 Charge ordering is 
indeed one of the fundamental routes to cause metal-insulator transitions in non-integer fill- 
ing systems. The nature of this transition, however, is not well understood yet. Here, we refer 
to non-integer filling as the state where the electron number per site is not an integer. The 
insulating phase of these systems usually have translational symmetry breaking in the charge 
degrees of freedom. This is in marked contrast to the case of integer-filling systems, which 
have integer numbers of electron(s) per site. Insulating phase in the latter systems, such as 
the Mott insulator phase, usually has no symmetry breaking in the charge degrees of freedom. 
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Here in this paper, we focus on systems at simple rational fillings that stabilize commensurate 
charge orders. We also focus on the interplay of the on-site and intersite Coulomb interactions 
in such charge-ordering transitions. 

There are several properties in charge-ordering transitions that require theoretical clarifi- 
cations: 

1) The system can become insulating with charge order but without the magnetic order: This 
is expected if one considers the strong-coupling limit. Systems with repulsive interaction at 
commensurate filling in this limit takes the lowest energy when the charge forms a certain 
spatial pattern to reduce the Coulomb interaction energy. In such phases, any excitations in 
the charge degrees of freedom has a finite gap. The system is then an insulator irrespective 
of the magnetic ordering, which may not be formed due to the quantum and/or thermal 
fluctuations. This expectation remains to be clarified except that limit. Experimentally, non- 
magnetic insulating phase is indeed observed at intermediate temperatures in wide range of 
real materials. 5 ' 6 

2) Many systems undergo metal-insulator transitions simultaneously at the point where the 
long-ranged charge order is formed i.e. no intermediate phase with charge order and metallic 
conductivity appears: Here, we do not in principle exclude the possibility of the charge-ordered 
metallic state. However, the above-mentioned feature is experimentally observed in most of 
the real charge-ordering materials. 3 ~ 5 ' 12-14 

On the other hand, using simple theoretical methods such as the mean-field approximation, 
one easily obtain charge-ordered metallic phase in a wide region of the parameter space as we 
discuss below. Mean-field approximations fail in reproducing the general tendency in realistic 
parameter region. 

Studies on metal-insulator transition induced by the strong electron-electron correlations 
require careful treatment of the two intriguing aspects of quantum mechanics, itineracy and 
locality. This is one of the most fundamental problems in condensed matter physics. In integer- 
filling systems, such transition is called the Mott transition 15 and has been the subject of 
intensive research. 7 These studies were based on a lattice model with only the on-site electron 
repulsion, called the Hubbard model. The Hubbard model is defined as 



Here, Ci a (c icr ) is the annihilation (creation) operator of an electron at an atomic site i with 
a spin index a. The number operator on site i with spin a is represented by rii a , while tij 
is an electron transfer from an atomic site i to j, and U is the local Coulomb repulsion. 
Recent development of numerical calculations has allowed detailed analyses on the Mott tran- 




(1) 



sition. 6 8 The dynamical mean- field theory and its extension have also applied to clarify 
finite-temperature properties of the Mott transition. 20 ' 21 
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On the other hand, as for metal-insulator transitions in non-integer filling systems, theo- 
retical mehods are still not well established. In this paper, we will study on the methods to 
treat those transitions induced by the charge ordering. 

Charge ordering in a broad sense includes the Wigner crystallization. 22 The charge order- 
ing in strongly correlated systems, however, is much different from the Wigner crystallization 
of the electron gas in the continuum space: First, there is a periodic potential formed by the 
positively (or negatively, if carriers are holes) charged ions, which may allow the ordering 
at relatively higher carrier concentrations. Second, as a consequence of the first, the order- 
ing occurs mostly at commensurate fillings and easily melts away from it. The importance of 
commensurability has been pointed out in ref. 23 In this paper, we call only this commensurate 
case as the charge order. In order to obtain a model for this charge order, we assume that the 
long-ranged part of the Coulomb potential is screened by the large carrier concentration. We 
define the simplified model, the extended Hubbard model as 



Here, Vij is the off-site Coulomb interaction. Other notation is the same as eq.(l). For sim- 
plicity, we restrict the range of V to only the nearest neighbor sites. We consider square lattice 
and hopping integral —Uj takes the value —t for only the nearest neighbor pair and zero 
otherwise. 

Charge ordering has been studied by several theoretical methods using the extended Hub- 
bard model defined in eq. (2). Here, we will examine some of the existing theoretical methods. 
The Hartree-Fock approximation (HFA) is the most standard method to treat the charge or- 
dering problems. 10 Using HFA, one can reproduce charge-ordered phase. On the other hand, 
HFA cannot reproduce the charge-ordered insulating phase without magnetic ordering. HFA 
also overestimates the ordered phases. The drawbacks arise from the neglect of dynamical 
and spatial correlation effects. Perturbative approaches 24 may be powerful tool for analyz- 
ing metallic phase of these systems. However, most of them cannot describe insulating or 
symmetry-broken phases. The dynamical mean-field theory (DMFT) can also be applied with 
appropriate mapping to the infinite dimensional model. It, however, fails in reproducing the 
insulating phase in the absence of the magnetic order. 25 This is ascribed to the fact that it 
treats the intersite interactions only in the static mean-field level. Exact diagonalization (ED) 
can also be a standard approach to this problem. 26 There is, however, finite size effect that 
becomes severe in determining metal-insulator transitions and/or charge ordering transitions. 
It is also difficult to obtain finite-temperature properties. 

Above considerations impose the requirements for an improved theoretical approach; (1) 
a non-perturbative method to describe the insulating and/or charge-ordered phases, (2) an 
analytic method that can describe the finite-temperature phases, (3) a method which can 
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equally treat both intersite and onsite interactions to realize insulating phases. 

In this paper, we analyze the charge-ordering phenomena by using the correlator projection 
method (CPM). 21 We show that it succeeds in reproducing the presence of the insulating 
charge-ordered phase in the absence of the magnetic order. 

The formulation of the correlator projection method is given in Sec. 2. Here, we summarize 
the properties of this method. CPM is a newly developed method that combines the equation- 
of- motion approach 27-29 and the dynamical mean- field theory. 19 It can treat the correlation 
effects beyond the mean-field level in a self-consistent manner. It is also a non-perturbative 
method that can describe the symmetry broken state and/or the insulating state. It consists 
of two parts: First, one derives the analytic form of the Green's function by the equation- 
of-motion approach. Second, one calculates the self-energy part by extending the dynamical 
mean-field theory in order to restore the fluctuation effects. These two parts are performed 
in a self-consistent manner on the basis of the operator projection theory and the dynamical 
mean-field theory. 

The dynamical mean-field theory is suited for treating the effects of dynamical quantum 
fluctuations. However, it may not be adequate to treat that of spatial fluctuations. CPM 
allows to overcome this drawback by taking account of spatial fluctuations. We concentrate 
on the system at commensurate filling in strong coupling regime, where the effect of spatial 
fluctuation is relatively small. 

The results obtained by the present approach are given in Sec. 3. CPM have reproduced 
charge-ordered insulating phase without order in the spin-degree of freedom. 

The result is in agreement with several experimental observations. Comparisons to real 
materials are given in Sec. 4. 

Sec. 5 is devoted to the summary of this paper. 

2. Formulation 

2.1 Correlator Projection Method 

As we have mentioned in the introduction, the correlator projection method (CPM) com- 
bines two procedures. 

The first procedure is based on the equation-of-motion decoupling approach originally 
named the operator projection method (OPM). 27-29 This approach, which was originally de- 
veloped for analyses of stochastic phenomena, was first applied to the Hubbard model by 
Roth. 30 She solved the equation of motion up to the second order and obtained a Green's 
function with two poles. This approximation is called the two-pole approximation (TPA). 
TPA has been further applied for the Hubbard model and some other models with strong 
electron-electron correlations. 31 ' 32 

Equation-of-motion approaches are suited for strongly correlated electron systems since 
one can treat the local constraint (the Pauli principle) correctly and can distinguish between 
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the double occupancy and single occupancy. Thus one can treat the strong on-site and off-site 
Coulomb interaction. The obtained Green's function correctly reproduces the high-energy part 
of physics. On the other hand, however, they tend to fail in describing the low-energy part of 
physics since they decouple the higher-order operators, which describe low-energy collective 
motions. 

The correlator projection method overcomes this drawback by introducing a refined self- 
energy at the point one truncates the equation of motion. When one solves the equation of 
motion for the Green's function, it leaves the self-energy as an unknown quantity. The equation 
of motion for the self-energy (, which we call the first-order self-energy T,^ from now on) is 
then given and leaves the second-order self-energy E( 2 ) as an unknown object. This sequence 
of continued fraction expansion continues. Using OPM, one can obtain the explicit form of the 
self-energy, which we call 'the higher-order self-energy', Y,( n \ The second procedure of CPM 
is to evaluate the highest-order self-energy, say Y,( n \ in self-consistent manner by extending 
the dynamical mean-field theory. The continued fraction is truncated at S^. Then, £( n ) is 
self-consistently determined by assuming that T,^ is momentum independent. 21 This reduces 
to the original dynamical mean-field theory when one takes n = 1. By solving a higher- 
order self-energy self-consistently, it systematically incorporates the spatial fluctuations with 
increasing n. For practical use, even n = 2 gives a significant improvement, where recovers 
momentum dependence. 21 

The first attempt of applying this method was given in the half-filled Hubbard model 
with the second-neighbor hopping i'. 21 It reproduced the metal- insulator transition at the 
parameters in good agreement with a reliable numerical calculation. 17 

2.2 Formulation of CPM 

Now, we summarize the concrete formulation of CPM. 21 

2.2.1 The first procedure; OPM 

One starts with a set of operators {4>i} for physical quantities. The equation of motion for 
the operators is expressed as 

—fr = ( 3 ) 

Through the remaining part of this paper, we take imaginary time representation; </>(r) = 
e rW 0e~ TW with the Hamiltonian of the system H. Here, we introduced operator u, which is 
defined from d>A = [A, H] for any operator A. The right hand side of eq. (3) is expressed as 

^ = E^ + #. (4) 

3 

The first term of the right hand side of eq. (4) is regarded as "system part" , but in practice 
describes the part within the subspace of {</>«}• Its coefficient is represented as e^ 1 ' 1 ), while the 
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explicit form is given later in eq. (5). The second term of eq.(4) describes the "environment 
part" which in fact describes the part orthogonal to the operator set {fa}. The decomposition 
into the two parts in eq. (4) is explicitly obtained by projection. We define the projection 
operator to the subspace of the operator set {(pi} as P\. It is defined for Fermionic/Bosonic 
operators as 



V / ii 



i,j v ' 13 

Here, cf) is the vector notation of {4>i} whose i'th element is fa. Then the system part is 
explicitly obtained from 

Ei"' 1 ^ = (5) 
3 

(2) 

The next step is to derive the equation of motion for <fi\ ; 

-^=^f. (6) 
By applying another projection operator, the right hand side becomes 

^ 2) = EW' 1) ^ + ^ 2 M 2) ) + ^ ) - ( 7 ) 

3 

The higher order equations of motion uj<j)f^ = Ya=i Ylj ^i™'^ <t>^ +<f>^ n+1 ^ are obtained similarly 
with the matrices defined as 

3 

Here, the n-th order operators are defined as 

P n x = {x, 4>^} ± ({{^ n \ 0(")t }±)) -i w 5 (9) 

n 

0("+l)= J](l-£)£>0(»). (10) 
1=1 

These equations generate (imaginary) frequency-dependent correlation functions. We in- 
troduce the notation ((• • -))icj n as the Fourier transformation of the correlation function 

«X;Yt }) ^ f dre^i-Xir)^). 
Jo 

We also introduce the Fourier transformation of the operators {fa} as 

i 

as well as the Fourier transformation of matrices A-ij as 



i 3 

where N is the number of sites in the system, and Vi is the coordinate of the i's site. The 
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equation of motion for the correlation function becomes 

tUniiPk ,9k lliun - e k +e k \\9k '^fe //«•>« +^fe \\9k '9k //«"»> I 11 ) 

and thus the Green's function becomes 
The self-energy £( n ) is given as 

4 n) = ((^ n+1) ;4 n+1)t )W- 

If one stops the expansion at the first step, and neglects T,^\ it obviously reproduces 
the Hartree-Fock Approximation (HFA). If one stops at the second step, and neglects £( 2 ), it 
gives the two-pole approximation (TPA) . 30 In case of the Hubbard model, for example, one 
can obtain the upper and lower Hubbard band at this step. Here, in the CPM, we calculate 
£( 2 ) to reach a better approximation. 



2.2.2 The second procedure; DMFT 

By means of OPM in the previous section, we obtain the Green's function of the system 



as 



Gk{iUn) = — — 1 i (13) 



lLU n - e k 



k ■ ' (2,2) v( 2 )/- -> 



k ^k 

In order to formulate the local approximation scheme, we reformulate the equation of mo- 
tion in the previous section by the path integral. In the following equations, we introduce 



the Grassmann field </>( n ) ) corresponding to the (normalized) operators cf>^/Ve^ n ' n 1 ) 
(<^>( n )t j\J e (n,n-i)) _ The partition function is given by 

Z = J n^ j (r)^,(T)e- s , (14) 

i 

where the action has the form 



s = / *- E msij + + E ^ yfW* + E 4 

E Wfr + M,- 2) )*f + E (^^# + v^| 2 >) . (15) 



<;3 

r 



1,3 i 

Here, the effective action is derived so that its functional derivative makes correct equation 
of motion for cj) and 9^ . If one could solve the action for 9^ exactly, one would obtain , 
where G^ is the full propagator for the <// 2 ) field. By using G^\ eq.(15) is rewritten as 



1,3 1 
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+ JdrJ dr'J2^\r)(-G\f ~\r - r'))4>f\r'). (16) 

In general, however, one cannot rigorously solve the action for fields. Here, we introduce 
local but imaginary-time dependent Weiss field to approximate the action for the 4>^ 
fields without the last two terms in eq. (15). Then, using the action is rewritten as 

s = J dr £ Mdrki + ^)<t> 3 + £ {^\Rf^ + #° y[W*) 

+ JdrJ dr'^\r){-Q {l) - l (r-r'))^\r') 

i 

+ £ (^ (2) V^f } + 0f } V^Vf } ) • (17) 

After integrating out (p^ field, we obtain 

S = Jdrj dr'Y^Hr) (-Q^ 1 (r - r')) . >,(r') + £ ($ 2 > \/^>f > + $ > 

(18) 

Here, the 'Weiss field' <?(°) is given as 

4° )_1 ( T - r ') = ("^Ai - e S' 1) ) <5 ( T " r ') " V^ 2 ^ 1 ^ - rOv 7 ^. (19) 
With the Weiss field, we can now solve Y,( 2 \ As was shown in the last section, E^ 2 ) has the 
form 

S( 2 ) = ((^);^)t))^, (20) 

where the field 4>^ can be expressed with the original <f> field. Since we have the Weiss field 
for (f) field, we can solve the 'second-order self-energy' in eq.(20) with some proper method. In 
this paper, we adopt the iterative perturbation theory (IPT) as we discuss below. 

The whole iteration procedure is explicitly given below: 
[1] IPT solution to £( 2 ) 
Evaluate the self-energy £( 2 ) as 

s(2) = ^E7^)((4 3) ;4 3)t )U. (2D 

k € k 

The operator cf)^ is evaluated by Wick expansion. For cj)^ expressed in the form <j>^ = 
^2 T k, k', qA'+q^k'^k-q, we can expand ^ by IPT as 

s<2) 4^ 7^,4<W^M% W - (22> 

k,k',q e k 

For the initial condition in the first iteration, one may use the Green's function for the free 

particle for Q^°\ 

[2] self-consistency equation 
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The local self-consistency condition is given as 



^>-*^V-*w <23) 

[3] Dyson equation for Q^> 

The Dyson equation for the field (f)^ is given as 

g^-\iu n ) = G[2 _1 (to; n ) + S( 2 )(iu; n ). (24) 
[4] Equation for 'Weiss field' G k 

From the eq.(19), we can calculate the momentum dependent Weiss field for the field fa as 

9f-\i^) = iu n - - ^Fa«(^)V^F- (25) 

After Fourier transform from Q^\iu n ) to Q^\t), Q^\t) is substituted to eq.(22) and we 
repeat the procedure from [1]. The procedure [l]-[4] constitutes the iteration loop. The pro- 
cedure is repeated until we reach the convergence. After the convergence, the Green function 
is obtained as 

G k {iu n ) = m . (26) 



iu n - e" k - 



i Wn -ef 2) -£(2)(iu, n ) 



2.3 Application to the extended Hubbard Model 

We now apply this method to the extended Hubbard model defined by eq.(2). Here, we 
introduce two methods in analyzing the model. 



2.3.1 Standard approach 

We first apply the OPM procedure in a straightforward manner. Now we take Cj CT as the 
first set of the operators fa in eq. (3). 
The equation of motion eq. (3) reads, 

WC iff = ^ 4]'^ C 3<r + VW- (27) 
3 

The matrix element is given as 

4j'^ = ^ij + (®ia)8i,j ~ Xij, 

where <Ev = Uni^ a + YljVi,j n j an d X^j = Vij(pj j i j(T ) are direct and exchange fields with 
the notation Pij j(T = c\ a Cj a , ni = + n^. We also define the following notations; 5&i a = 
$ia ~ {<S>ia}, $Ya = Ylj V ij n j- $®Ya = ®Ya ~ i®Ya) ■ The second operator i{; ia is obtained as 

We can now proceed to the equation of motion for the second set of operators, which is 
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formally written as 

^ = EH- ,1) ^ +c S ,2) ^)+^ ( 28 ) 



where the matrices are written as follows 



5 m 



i 2,1) = m^j}) = {mf)k 3 - m Pi ,i)\\ 

$ 2) = -iij-Xij + ((**)-8fH)Sij, 
tij = (ut^SmSm + Si S t - Af Af) + {6*Ys*Y?) Ui 
x((^))-%, 

-^E - E ^-'^) I ( x - 2 < n -» ( e(2,1) ) _1 ' 

i .?>' / 

^ = V t3 (p 3ta )(25<5> t(T 5<$> 3a - {5<S> ia f - (5<S> 3a ) 2 }, 

{ ia) m la Y) ■ 

Here, Si is the spin operator on the site i, defined as Si = ^c\ a cr aj pCj/3 with Pauli matrices 
er. The singlet pair operator on the site i, Af is defined as Af = Cj-fC^. The motion of the 
interacting charges is reflected in 5m The term with $ cfr comes from the interaction part. It 
differs from ($) because of the Pauli principle, (e.g. ni^ a ni- a Ci a = rii- a Ci a ) 

The remaining operator, ( the third operator of the right hand side of the eq. (28) ) R 
reflects the fluctuation of <5<E>. It is given as 

R%u = U ^ ^ (( tij) (pij^—v Pj,i,—o) Cia tij5rii^ u Cj a ) 
3 

jcr> V I ) 



+ + 6p5 io)^j<T + R 'ic 



Here, R' icr contains higher order operators as 



Ricr — 5Qi a Ci a 

with 5Qi a defined as 5Qi a = (5&i a ) 2 — j^- — la \l 5&io- The correlation functions of this op- 
erator give a convolution of Green's function with charge/spin correlation functions of more 
than three terms such as {{5rii5nj5ni)) in which indices i, j and I are all different. Assuming 
that these higher order fluctuations are small, we neglect this operator R' ia . 
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Here, in the matrix elements, one finds equal-time correlation functions like ((5&) 2 ) that 
cannot be uniquely determined by the single-particle Green's function. There are several ways 
to determine these correlation functions. 

(1) The second-order perturbation theory: One way is to make use of one of various perturba- 
tive methods. Among them, we adopt the second order perturbation theory. Here, we do not 
start from the free-particle Green function, but use obtained (renormalized) Green function 
and take up to the second order diagrams to include correlation effects. Although it is one of 
the simplest approximations, it is favorable for our purpose here since it conserves the local 
summation rule; 33 

iEE(^w^f p) ( iw »)) = w(2-w). (30) 

u-n, q 

Here, Xq° h ' \iu n ) = ((n q ;n^ q )} iuJn and Xq P \^n) = (((%r ~ n gl); ( n -<?t ~ n -g|)))^n are tne 
charge and spin correlation functions, respectively. By the second order perturbation theory, 
they are derived as 

Xq ch) (iu; n ) = 2{ x ° q (iL0 n ) + x° q (i"n)(U + 2V q ) X ° q ^ n ) 
k,k' 

Xq SP) (iuJn) = 2{X° q (iu> n ) - X q (i"n)Ux q (iUn) 



1 ( 32 ) 

T7f 2^ Xfc; g («^n)Vfc_fc'Xfc'; g («^n)}- 



N 2 
k,k> 

Here, Xqi^n) an d Xk-qfan) ar e zeroth order susceptibility functions defined as 

X q (iu n ) = Jj-j^2Y^ °k+q( iuJ i + i^n)G k (iuJi) (33) 

k 

Xk; q (iu n ) = g ^2 G k+q {iui + iu n )G k (iuJi), (34) 

where Gk(ito n ) is determined from eq.(26). On the other hand, it may underestimate the local 
correlation effects especially in the case of strong coupling. 

(2) The exact diagonalization (ED) of small cluster: We can also take a completely different 
approach by simply replacing the correlation functions with the result of exact diagonalization 
(ED) of a small cluster. Adopting ED, one can take account of the strong local correlation 
effects. Here we emphasize that we only need short-ranged equal time correlation functions 
that are not likely to severely suffer from finite size effects. In order to treat the symmetry- 
broken state in the ED calculation, we divide the lattice into two sublattices, A and B, and 



11/28 



J. Phys. Soc. Jpn. 



Full Paper 



apply small symmetry-breaking field h as follows; 

~H c s = Wehm - h J n i ~ ^ n j 

\ieA jeB 



Here, the first term TCerm is the original Hamiltonian of the extended Hubbard model given 
in eq. (2), while the second and third term indicate the symmetry breaking field. The field 
h is set h = O.OOli. The shortcoming of this approach is that the correlation functions are 
determined independently of the lattice Green's function, and are not fully self-consistent. 

We adopt the former approach to derive the charge ordering transition where the con- 
sistency to the Green's function is particularly important. On the other hand, we adopt the 
latter approach to analyze the detailed properties for the charge-ordered phase. 

2.3.2 Formulation based on the CDW basis 

In order to further analyze the properties of the charge-ordered phase, we introduce a 
different approach using different basis representation. Although in principle, the operator 
projection theory does not depend on the choice of basis, in actual calculations, however, 
details of the results do depend on the chosen basis set due to the truncation of the operator 
expansion at a finite step. Because of the truncation, it is better to employ the basis function 
which is the eignfunction of the same symmetry as that of the correct phase. For example, 
in the charge-disordered state, the basis function of the free-particle Hamiltonian gives faster 
convergence and gives better results at a fixed level of the truncation. This is the case in the 
previous subsection, namely, the standard approach. On the other hand, if the system is in the 
charge ordered phase, the basis function of charge density wave (CDW) mean-field solution 
is expected to give better results. Here, we introduce an approximate formulation which is 
valid under large polarization. The idea is to solve the equation of motion for the operators 
represented by the eignfunctions of the mean-field Hamiltonian for the charge order. 

In the following equations, we consider the square lattice as an example. Hence, at quarter 
filling, two sublattices A and B are enough to discuss possible charge ordering. The dispersion 
in the two-dimensional square lattice is denoted as — tk = —2t(cos(k x ) + cos(k y )), and electron 
density per site as (ua) and (ne), for each sublattice. We also use the notation (n) = \{{ua) + 
(ub)) We take new basis as 



Here, c^(c fc T ) is the annihilation (creation) operator of the Bloch state at sublattice R (R = 
A or B), with momentum k . Here, and in the following equations in this section, k refers 
to the wavevector inside the reduced Broullin zone (RBZ). We fix Uk and Vk as the eigenvector 
of the mean-field Hamiltonian matrix e HF . The Hubbard term £/^i n iT n il is treated in later 



Ofc = U k Ck + V k Ck 

bk = -v k ci + u k Ck- 



(35b) 



(35a) 
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calculations. We define the CDW order parameter Acdw as 

Acdw = 2 x W((n A ) ~ (n B )). (36) 

The order parameter is to be determined in a self-consistent manner. The mean-field Hamil- 
tonian becomes 

e HF = ( iV ( n ) ~ A cdw ~t k \ 

\ -t fc 4V(n) + Acdw- / ' 

Hence Uk and Vk is obtained by the following equations: 

e H F f u k v k \ = ( u k v k \(\£ \ m 
V -Vk u k J V -v k u k J \ Af / 

Here, A£ = 4F(n) - A fc and Af = 4V(n) + A fc are eignvalues of the mean-field Hamiltonian 
with 

Afc = ^A CDW + ^ . (39) 
The eigenvectors are determined from 



Uk = (40a) 



'a 




+ Acdw 


2A fc 


'A 




- A C dw 



Vk = fa, • < 40b) 

We perform inverse Fourier-transformation of eq.(35) to obtain site representation as 

cti = — U V afce* kri i £ A, (41a) 

. / AT ^ — ' 



N k 

b i = ^E 6 ^ k rj ( 41b ) 
^ k 

The last expression also satisfies the anticommutation relation; 

{ai,a\,} = <5 M /, (42a) 

{6^,6],} = ^, (42b) 

and any other pairs than these anticommute. Note that in spite of the site index, they have 
nonzero amplitudes over many lattice sites as 

V j 

b j = -Yl + Yl u 3j'4 ' ( 43b ) 

i j' 

where Uij = u(ri — rj) and Vij = v(ri — rj) indicate the Fourier transformation of Uk 
and Vk- The Hamiltonian can be written with the new operators. We introduce the following 
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notations; 



Si.x, — 



Qj,x 



-v iiX x e B, 

Vj, x x £ A, 

Q>x,<j X G -A, 

b x ,a x e B. 



We write the Hamiltonian as 

H = H m{ + Hu + H'. (44) 

Here, H m { is the mean-field terms; 

H mi = - f iAa C ja + 22 V ij( n i) n i ~ \ 2 2 V i,j( n i)( n j) ( 45 ) 

* 3 i 3 

and K{/ is the Hubbard interaction terms; 

H U = Y J Unf^ + 2 Unf^. (46) 

The rest part of the mean field Hamiltonian is given by 

n' = J2 v i,J Sn f Sn f- ( 47 ) 

(i,3) 

Then, H m { is rewritten as 

«mf = 2 A M' a i a ^' + 2 A M 6 i' + ( 48 ) 

i,i'eA j,j'eB 

Here, A^ and \? m is the Fourier transform of A^ and A|? within each sublattice, respectively, 
and E m f is the mean field energy; 

A 2 

E mi = 2V(n)N e + -gpiV, 
with N e being the total electron number. Whereas H' and Hu are written by substituting 

T^ia — ^ ^ Si,xSi,yd xcr dya (49a) 
x,y 

n i° = 2 Qj,*Qj,v d L d w ( 49b ) 

into the original number operators n xa in eqs. (46), (47). 

In general case, the above-mentioned procedure yields complication. However, it becomes 
simple in the largely disproportionated case. In the following of this section, we consider such 
case. We assume the CDW order parameter Acdw to be considerably large so that we define 
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an expansion parameter p as 

P 2A CDW ^ ^ 

Then and v^, in eq.(40) are expanded in power of p up to p 2 as 

u k = l--J^ + 0(p 4 ) (51a) 
4 ^CDW 

v k = ^r^ + 0(p 3 ). (51b) 
^CDW 

Up to this order, their transformations to the real-space representation have nonzero values 
only within the nearest neighbor sites. We introduce the notation pij as 

i p (i,j) for the nearest neighbor pair 

Pi,j = \ 

[ otherwise . 

Then u and v become 

Ui,i> = 6 iti > +0(p 2 ), (52a) 
Vij=Pij + 0(p 3 ). (52b) 

In the above equation, i and %' are on the same sublattice, whereas j is on the other. The 
original operators cf and are expressed by Oj and bj as 

4 = fiai-^Pijbj), (53a) 



with / being the normalization factor / = , 1 The number operators are rewritten as 

V 1 + 4 P 



"f 



where 7?f CT = a| CT ai CT , and 77^ = bj a bj a . At this point it is clear that the factor / can be absorbed 
in the interaction parameter U and V as U = Uf 4 and V = V f A . The Hubbard interaction 
term contains 'off-site' terms that describe the interaction with states centered at different 
sites. 
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We write down the kinetic equation for the operator ai up to the first order of p as 
ua ia = ^2 \ya Va + A$ ia a ia - ^^(A^v - A$ ja )b ja 



j 



(55) 



V 

j j 

= ^ Ki' a i'o- + U(rj? a )a ia 
V 

j 

with the "second" operator fi a , defined as 

f ia = 5$i a a ia - ^2pij{5$ ia - 6$ ja )b ja . (56) 

j 

Here, we have introduced the notation A$j CT , defined as A§i a = <3?j CT — Ylj^ij( n j)- Other 
notations are the same as in the previous subsection. In the derivation of eq. (55), we have 
used the fact that terms which are off-diagonal in 'a' and '6' bands such as (a\ a bj) are of order 
p or higher, and terms like p{a\ a bj) can be neglected when we consider the order up to 0(p). 
The equation for the fields bj a 's is also obtained in the same way. 

^ b i° = ^2 X i,™ b ™ + A$ ja b ja - ^2pji(A^ ia - A$ ja )a ia 

m i 

= Yl X 3,mbma + U {v)-a)bja (57) 
m 

i 

with the corresponding "second" operator gj a , 

9j<r = 5<& ja b ja - ^2pji(5<£ ia - 5<£ ja )a ia . (58) 
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The equation of motion for the field fi a becomes 

3 3 

= Cb{<&i a )a ia + 5<S> ia (Y^ + A<Z> ia a ia ) 

n 

~ ^Piji^ia ~ <S>j*))bja 

3 

~ J2pij(5<$> i<T - 5<$> ja )(^2 X jrnbma + ^jabja) 
j m 

_ ST^ >a)(2,l) >b)(2,l), 

~ 2-^i ™ n<T ^ Z_> ij 3° 

n j 

sr^ >a)(2,2) f (a6)(2,2) j 

n j 

The equation for Qjcj is also written in the same way; 

r , n . (6a)(2,l)_ V^,(»)(2,1)l 

n m 

, V" f (H(2,2) f , ,(»)(2,2) o 

j m 

The explicit forms of the remaining operators and will be given later. The matrix 
e ( 1 > 1 ) i s defined as 

^(1,1) = U {v a_ a)5in + XL (61) 

e (f >(M> = -py(Lty£_ ff > - C/(^_ CT )), (62) 

whereas e^ 6 ^ 1 ' 1 ) is obtained from g( aa )( 1 ' 1 ) by replacing 'a (A)' with 'b(B)' and p with —p. 
Here, e^"^ 1,1 -* is obtained from the symmetry property e^"^ 1 ' 1 ^ = e^ 6 ^ 1 ' 1 ^. The matrix e 1 - 2 ' 1 ) 
is also derived up to the first order of p as 

£0(2,D = m) 2 )6 ^ (63a) 
e (a 6 )(2,l) = _ ( ^. (( ^ a) 2 _ (5^)2)). (63b) 

Here, however, the matrix e*- 2 ' 1 ) need to have positive eignvalue since it is the norm of the 
operators fi a and gj a . We modify it within the second order of p to guarantee the positivity. 
It is modified, in the momentum representation, as 

(2,d i ( i Pk\( ((s*i*) 2 ) o W i - Pk \ 

€k 1 + <A)\-PH 1 ){ <(M>, CT ) 2 > ){ Pk 1 J' [b) 
where p^ is the Fourier transformation of pij . Since the right hand side of eq. (64) is a unitary 
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transformation of a diagonal matrix with positive matrix elements, e^ 2 ' 1 ) in this form has 
non- negative eignvalue. Note that eq.(64) is equivalent to the previous eq.(63) up to the first 
order of p. In order to write down e^ 2,2 ) , we define the matrix vti^y^ , which is related to the 
matrix e^ 2 ' 2 ) with 

e (2,2) =m (2,2) (e (2,l) ) -l. (65) 

The matrix m^ 2 ' 2 ^ is given by 

mt a){2 ' 2) = {u(l- (V-.)) + A^ ) e !r )(2 ' 1) 

+ {{5$> ia f)8 in . (66) 

m (a 6 )(2,2) = _ pij{{% _ _ p .. {T . § ^. a _ 

i 

n 

{A<f> ia 5 in + \f n )8<S> ia ) (67) 
In these equations, we have introduced the following quantities; 

1? = UJl a + %{Jl + JJ_ a ), (68a) 

3 xa — ^xyi^xadya ~ dy a d xa ) 

"F Pxy ( a xcr^yzbzv ~ ^yz^la a xu 

+ b\ a \x Z a Z(T - \ xz a\ a b y(T ), (68b) 
5fi r = (T xa 6<f> xa {l - 2n la )) 

= (\J 2 K x -„ + VZyVCy* + Ky- a )) (1 - 2{n ia )), (68c) 

Kc = E Xr *y(4*d y „), (68d) 

y-.yj^x 

C r xya ee (8*%6tf) + U 2 {\5rfjrf y + S x S y - A*Ag>. (68e) 

Here, the label 'r(f)' indicates either l A(By or l B(A)\ whereas =F takes — sign if the index 
x in the left hand side is in the l A' lattice, and + if x is in the l B' lattice. Greek indices 
fi, v represent both site and spin indices. From eq. (65) is also obtained by replacing 

a(A) with b(B) and p with —p. The matrix 

m (ba) ig 

obtained by the symmetry relationship 
m^f 1 = mlj b \ By using these quantities, the remaining operator Lj CT in eq.(59) is expressed 
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as 

L ia = T ia a ia + 5$ ia ^ >Hn a na ~ ^jPij{% - Tj)b ja + L \ a - ^ ^tnfna, (69) 
n j n 

where X A is defined as X A n = ^/ {-8n A 6u +C A \ A ) ((e^ ' ^ 2 ' 1 )) -1 ^ . Here, L' ia comes from 
higher order fluctuation terms as 

L' ia = 5Q ia a ia - ^2pij(SQ ia - SQ ja )b ja . (70) 

j 

Here, SQ is defined as 5Qi a = (<5$j CT ) 2 — wf^^L? <5&ia as the same as in the previous subsec- 
tion. We ignore them in the same manner as in the previous subsection. The operator Rj a is 
also written in the same way as 

Rj<r = Tjabja + 8®ja ^ Xf m b ma - ^^Piji^i ~ %) a iv + R-'ja ~ ^ ^fm9ma- (71) 

mi m 

Here, is given as 

R'ja = SQjabja ~ ^Pji^Qia ~ 5Q ja )a ia , (72) 

i 

and it is neglected here. 

This formulation is valid only if the CDW order parameter is large enough to ignore the 
higher order of the expansion parameter p = 2 Ac DW • However, we expect that this formulation 
gives better result even when the order parameter amplitude is modest since the CDW basis 
becomes correct in the strong coupling limit up to the second order in y. If the intersite 

interaction is large so that t <C Aqdw, one can expand the quasiparticle dispersion eq.(39) 

t 



up to the second order in -i as 



t 2 



Xk — A-CDW + T • 

Performing Fourier transformation, one can obtain the effective hopping as 

where faj = -jj e lk ^ ri ~ r ^ (2(cos(k x ) + cos(k y )) 2 is the next-nearest neighbor hopping. This 
corresponds to the second order perturbation around the charge-ordered state. We refer to 
this formulation as 'CDW basis approach' in contrast to the formulation introduced in the 
previous subsection 'standard approach', which takes {c^} as the operator basis. 

3. Numerical Results 



In this section, we show the result of the numerical calculations. 

The following calculations have been performed on the 32 x 32 lattice with Matsubara fre- 
quencies 2048. The spectral functions have been obtained by analitic continuation using Pade 
approximation. 34 The obtained real frequency data are smoothened by taking convolution 
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Fig. 1. Charge and/or spin distribution of the ordered patterns obtained by CPM (with CDW basis). 
Large circles indicate charge rich sites, while arrows indicate spin polarized sites. Left panel shows 
the solution for T < Tn with magnetic order. This pattern is the same as the HFA solution (3). 
Right panel shows the solution for T > T/v without magnetic order. 



with the Gaussian distribution function N(lo: a) = ,— e ^ with a = 0.05. 

Through this paper the energy unit is set as t = 1. First, we show the phase diagram of the 
extended Hubbard model obtained by CPM. We first derived the charge-ordering transition 
line by the 'standard approach', taking Cj as the operator basis. We then have changed the basis 
to the 'CDW basis' in order to further analyze the region where the CDW order parameter 
Acdw is large. 

First, we compare the phase diagram obtained by CPM with the corresponding one ob- 
tained by the Hartree-Fock approximation (HFA). Here, the both calculations were performed 
allowing 2x2 sublattices. In the HFA calculation, we have examined longer-ranged ordering 
patterns allowing up to 4 x 4 sublattices but no other type of order was obtained as a stable 
phase. One clear difference of the present CPM result from the corresponding result by HFA 
is seen in the existence of the insulating phase with charge order but without the magnetic 
order (non-magnetic insulating phase). Here we note that, even in the CPM calculation, an- 
tiferromagnetic order is indeed obtained at low temperatures. The obtained ordering pattern 
is shown in Fig. 1, which is in agreement with strong-coupling analysis as well as the phase 
(3) of the HFA phase diagram. At first thought, by taking account the fact we are considering 
two-dimensional system, we can argue that this result violates Mermin- Wagner theorem. 35 
This violation is, however, attributed to the mean-field like nature of the calculation method. 
We can expect that the problem would be reduced if we take larger sublattices. 

Here, we emphasize that in the CPM calculation, magnetic order is suppressed to low 
temperatures. In Fig 3 we have shown the Neel temperature for several parameters. This low 
ordering temperatures indicate that the magnetic order is induced by a small energy scale 
J e ff, which is the exchange energy in the charge-ordered phase. The explicit calculation of this 
energy will be given in the Appendix B. Here, the calculated Tn is almost the same order as 
the exchange energy scale obtained by the strong coupling expansion. Although Tn is higher 
than J e ff especially in the weak coupling region, this is atttributed to the deviation from strong 
coupling assumption. On the other hand, metal-insulator transition as well as charge-ordering 
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2 4 6 8 10 12 14 



Fig. 2. Phase diagram for the extended Hubbard model at quarter filling. The left panel shows 
the result obtained by CPM, while the right panel shows the corresponding result obtained by 
HFA. Temperature is set T = 0.15. The dotted lines indicate continuous transitions, while the 
solid lines indicate discontinuous transitions. In the left panel, COI, COM, and CDM indicate the 
charge-ordered insulating, the charge-ordered metallic, and the charge-disordered metallic phases, 
respectively. In the right panel, the labels (0) to (4) indicate different ordering patterns of spin 
and/or charge obtained by the HFA calculation. Each ordering pattern is shown in Appendix B. 




10 11 12 13 14 
U 

Fig. 3. Plot of the Neel temperature Tm as a function of interaction parameter U, V. The result is 
obtained by CPM based on CDW basis. 



transition occur in larger energy scale. We have confirmed that the charge order as well as the 
metal- insulator transition up to T m 0.5, which is several factor higher than Tjv- This is in 
contrast to HFA, in which magnetic long-range order is necessary to reduce double occupancy. 
Thus magnetic order is induced by larger energy scale U, and it exists even in relatively high 
temperature. 

Using CPM (CDW basis approach) and taking account of short-range correlation effects, 
we can obtain a phase with small double occupancy but without magnetic long-range order. 
This correspond to the non-magnetic insulating phase. Using HFA, the square lattice system 
at quarter filling become insulating only if there is a spin symmetry breaking which splits the 
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Fig. 4. Local density of states p(u>) of the extended Hubbard model at quarter filling. Temperature 
is set T — 0.15, while the interaction parameters are set V — 2.5, U = 4.0, 12.0, and V = 3.5,{7 = 
4.0, 12.0 respectively. Solid line indicates the result obtained by the CPM based on standard basis, 
while the dotted line indicates the result obtained by the CPM based on CDW basis. 



band into four. On the other hand, by taking account of dynamical correlation effects as in 
CPM, one can reproduce insulating phase without magnetic symmetry breaking. 

We note that application of CPM in the 'standard approach' does not reproduce the metal- 
insulator transition. This is because the operator projection analysis up to the second-order 
expansion cannot fully reproduce the correct dynamics around the charge-ordering transition. 
On the other hand, the formulation based on the CDW basis reproduces the metal-insulator 
transition. This is attributed to the fact that this formulation recovers the strong coupling 
limit up to the second order in y as we have seen before. This basis function gives qualitatively 
correct results even at the truncation up to the second order. 

Here, one should note that the two approaches do not contradict at intermediate cou- 
pling since they are connected with (approximate) unitary transformation. Therefore, one can 
interpolate these two approaches and can obtain qualitatively correct results in the whole 
parameter space. In Fig. 4, we have plotted the spectral function obtained by the both ap- 
proaches. The two spectra show similar structure for intermediate amplitude of V. 

In Fig. 5, we show the spectral function of the metallic as well as the insulating phase. 
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Fig. 5. Local density of states p(u) at the metal-insulator transition point, The result is obtained 
by CPM based on CDW basis. Left panel shows the result for T = 0.15. The parameter is set 
as V = 2.5 T — 0.15, and U — 6.0, U — 6.1 respectively. The solid line indicates the metallic 
solution while the dotted line indicates the insulating solution. The right panel shows the result 
for T = 0.50. 
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Fig. 6. Order parameter of the charge-ordered phase (An) plotted as a function of on-site interaction 
[/.Temperature is set T = 0.15 (left panel) and T = 0.50 (right panel), respectively. 



The former phase has a large density of states at the Fermi level, while the latter phase has 
small density of states, which becomes vanishing in the zero temperature limit. 

Our result shows that the transition is of the first order, which is typical to (finite temper- 
ature) metal-insulator transitions. The order parameter of the CDW state (An) = (ua) — («b) 
jumps at the transition point as is shown Fig. 6 

At fixed V, the order parameter (An) decreases as a function of U in the metallic phase, 
but it takes nearly a constant value in the insulating phase. This is explained as follows; in the 
metallic phase with delocalized wave function, charge disproportionation increase the double 
occupancy and thus is disfavored in terms of the onsite repulsion term. On the other hand, in 
the insulating phase where the charges are localized each per one site, charge disproportion- 
ation does not severely increase the double occupancy. Moreover, suppression of the charge 
fluctuations results in larger charge disproportionation in the insulating phase. 
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In this picture, however, the insulating phase occur only in the charge-ordered phase. This 
means that there exist charge-ordered metallic phase for a certain parameter region. On the 
other hand, such phases are not widely observed in experiments. This issue will be discussed 
in the next section. 

4. Comparison to experimental observations 

In this section, we compare the obtained results with the experimental observations. 

As we have mentioned in Sec.l, charge ordering occurs in various types of strongly cor- 
related electron systems. Among them, charge ordering induced by the intersite Coulomb 
interactions is considered to be realized in a certain class of organic conductors. 10,11 Here, 
we note that there are number of transition metal oxides (TMO) that shows charge ordering. 
There are, however, some difficulties in comparing those materials with this result since in 
many cases, TMO's have more than one band and the spin degree of freedom is strongly 
coupled to the orbital degree of freedom. The result obtained here, by simple one band model 
with only electron-electron repulsion, may not be applicable. On the other hand, many of 
organic conductors have single conduction band and the result of single-band model is di- 
rectly applicable. In addition to this, these materials provide good examples where the charge 
ordering temperature is much higher than the temperature at which spin degree of freedom 
is lost (magnetic-ordering temperature, spin Peierls temperature etc.) as we will see below. 

For example, #-(BEDT-TTF) 2 MM'(SCN)4(M=Rb,Cs, M'=Zn,Co) have large two- 
dimensional anisotropy and strong electron-electron correlations compared to the band- 
width. 10 ' 11,36 ' 37 In addition to this, they have weak dimerization and regarded as quarter 
(hole-) filling system. Thus they offer good candidates of the charge-ordering system. In fact, 
6>-(BEDT-TTF) 2 RbZn(SCN) 4 , for example, undergoes MIT at T M n = 195K at ambient 
pressure, where the electrical resistivity shows sharp increase, with no sharp change in the 
magnetic susceptibility until a magnetic transition to spin Peierls state at Tgp = 50K . 5 The 
evidence of charge disproportionation was reported from the data of 13 C NMR shift and 
relaxation time measurement. 14 Between 195K and 50K, an insulating phase without any 
magnetic order is observed, which is in agreement with our numerical calculation. The dy- 
namical as well as spatial correlations are considered to be important in the above-mentioned 
system. The remaining magnetic degrees of freedom is lost through the spin-Peierls transition 
in the low temperature phase. This transition is considered to be induced by the coupling to 
the lattice degrees of freedom. Here, however, we do not further discuss this problem since 
we do not take account of lattice degrees of freedom. In this paper, we have adopted purely 
electronic model in order to study the physics in charge degrees of freedom. Study of the 
electron-lattice interplay remains as a future work. 

Here, it should be noted that the ordering pattern observed in 6— (BEDT- 
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TTF)2RbZn(SCN)4 is not the checkerboard order but of the stripe order. This difference, 
however, may be attributed to the difference of the lattice shape. In the real system, there 
exist finite values of the Coulomb interactions as well as the hopping in the diagonal directions. 

More important difference is the existence of the charge-ordered metallic phase in the 
numerical calculations, while it is not widely observed in experiments. This issue has been 
pointed out in Sec.l. One possible reason for the difference is the overestimate of the charge- 
ordered phase in the numerical calculation. Although CPM takes into account spatial short- 
ranged correlation effects, it still has a mean-field like character. This may result in the 
underestimate of spatial fluctuations and overestimate of charge order. If one could treat the 
spatial fluctuations more accurately, the phase boundary between the charge-ordered phase 
and charge-disordered phase may merge to the phase boundary between the metallic and 
insulating phase. On the other hand, however, this problem might also be attributed to the 
limitation of the extended Hubbard model itself as a realistic model of real materials. In 
real materials, charge ordering often accompanies some additional symmetry breaking such 
as lattice distortion or dimerization, which may further stabilize the insulating phase. 

5. Summary 

We have investigated the two-dimensional extended Hubbard model with the correlator 
projection method (CPM). CPM is a newly developed method which can treat correlation 
effects in a self-consistent manner. 

Applying CPM, we have shown that a metal-insulator transition is induced by charge 
ordering. The insulating phase has antiferromagnetic order below T/v, which is as the same 
order as the exchange interaction energy J e g , but is non-magnetic insulator above Tjy. This is 
consistent with the strong-coupling picture as well as several experimental observations. The 
temperature scale of T/v is several factor smaller than the charge-ordering temperature Tco 
and/or metal- insulator transition temperature Tmit- We also emphasize that the obtained 
result is in contrast to the Hartree-Fock approximation. The difference of energy scale as 
well as the existence of non-magnetic insulating phase is obtained by taking account of the 
dynamical and spatial correlation effects. 

However, we have also obtained a charge-ordered metallic phase that is not widely observed 
in real materials. The difference may arise from the insufficient treatment of spatial fluctuation 
effects. More accurate estimate of short-ranged spatial fluctuation effects is left for future 
studies. The emergence of charge-ordered metal in the present result in general contradiction 
to the experimental observations may also be attributed to the limitation of the extended 
Hubbard model itself, where the coupling to lattice is ignored. 

Although we clearly need further improvement, we have succeeded in obtaining a new 
formulation that can treat metal-insulator transitions induced by the charge ordering. The 
result also indicates that the qualitative aspect of the experimentally-observed metal-insulator 
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transitions in charge-ordering systems can be reproduced by a simple model as the extended 
Hubbard model if one carefully treat the correlation effects. 
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Appendix A: Phases obtained by the Hartree-Fock approximation 

The pattern of the charge and/or spin obtained by the Hartree-Fock approximation is 
shown below. We have examined possible spatial pattern up to 4 x 4 period and obtained the 
following 5 phases as stable ones. We do not exclude the possibility of obtaining lower energy 
state by allowing spatial patterns with longer period than 4x4. 
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Fig. A4. Charge and/or spin distribution of the ordered patterns obtained by HFA. Large circles 
indicate charge rich sites, while arrows indicate spin polarized sites. (0): no order (l):(ir, 0) spin 
order (2):(ir,n) charge order (3):(7r, it) charge order with (n,0) (0, 7r) spin order (4) ferromagnetic 
order 



Appendix B: Derivation of exchange coupling parameter 

The exchange interaction between localized charges in the checkerboard charge ordered 
state is derived from virtual processes shown in Fig. Bl. We define the interaction strength 
parameter P = mm(U, V). We first consider the strong coupling limit where P/t is infinitely 
large. Then we take account of the effect of finite value of hopping t by order by order in t/P. 
In Fig. Bl, the second order process (b) merely adds a constant energy shift irrespective of 
the spin configuration. Non-trivial effect comes from the fourth order perturbation. Contribu- 
tions from this process differ depending on spin configurations. In order to consider effective 
exchange interaction in the charge ordered phase, we introduce y/2 x y/2 magnetic lattice. The 
effective interaction is then expressed in terms of the magnetic lattice as J c g in the nearest 
neighbor (n.n.) pairs and J' cS in the next-nearest neighbor (n.n.n.) or the diagonal directions. 
Each is obtained as, 

16t 4 1 1 J_ 

J cff ~ ^72 \Jj + 4V + U + 4^)' 

4/ 4 1 2 
eff 9V 2{ U AV + U' 
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Fig. Bl. Virtual hopping processes up to the fourth order in t/P. Circles indicate charges with 

no spin specification, while arrows indicate spin.(a)the ground state configuration in the strong 

coupling limit (zeroth order in t/P). (b) the second order perturbation process that gives uniform 
t 2 

energy shift 4—. (c) spin pair of the second-nearest neighbor sites and (c'l) to (c'3) indicate 

different processes that contribute in the the fourth order perturbation. The contribution from 

At 4 4i 4 4i 4 

each process is ^v) 2 (W + U) ' (3V) 2 U ' &nd (3V) 2 (4V) ' respecivel y- ( d ) thc spin pair of thc 
third-neighbor sites and (d'l) to (d'3) indicate different processes that contribute in the the fourth 



order perturbation. The contribution from each process is 



4t 4 



and 



respecively. (e) The effective y/2 x y/2 magnetic sublattices for the charge-ordered phase. The 
effective exchange is described as J e g and J' cS . 



respectively. Here, the effective Hamiltonian has the form 

T~L = Jcfi Si'Sj + J'cS ^ Si-Sj, 

<i,j>Gn.n. <j,j>6n.n.n. 

where S is the spin- 1/2 operator at the charge-rich site of the charge-ordered phase, and 
Yl<ij>enn indicates the summation over all nearest neighbor pairs, while J2<ij><=nnn m " 
dicates that of all next-nearest neighbor pairs (in the sense of y/2 x y/2 magnetic lattice). 
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